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The  Equation  of  State  for  Anunonia* 
by 

Lester  Haar  and  John  Gallagher 

1.  Introduction 

This  report  has  been  prepared  in  answer  to  several  urgent  requests  for 
detailed  computer  programs  for  the  thermodynamic  properties  of  ammonia.  In 
it  we  present  a  description  of  an  extensive  correlation  of  the  thermodynamic 
properties  of  ammonia  published  elsewhere'^,  along  with  basic  equations  for 
the  description  of  the  properties  of  ammonia  and  detailed  listings  of 
computer  programs  based  on  these  equations.    With  this  information, 
thermodynamic  properties  can  be  calculated  for  temperatures  ranging  from 
the  triple  point  temperature  to  about  5/3  the  critical  temperature  and  for 
pressures  ranging  from  the  dilute  gas  to  about  8000  bars.     The  reference 
state  for  all  properties  is  the  ideal  gas  at  zero  kelvin.     The  physical 
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constants  used  are  consistent  with  those  recommended  by  Cohen  and  Taylor  . 
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The  mass  of  a  mole  of  ammonia  was  taken  to  be  17.0306  grams  . 

2.  The  derivation  of  the  thermodynamic  surface. 

In  the  correlation,  selected  P,p,T  data  were  fitted  to  an  analytic 
equation  of  state  in  the  least  squares  sense.    A  detailed  discussion  of  the 
process  of  data  selection  and  of  selection  of  parts  of  the  analytic  equation 
is  contained  in  reference  1.     The  equation  so  obtained  contains  the  pressure 
as  a  44  term  double  power  series  function  of  temperature  and  density. 
The  resulting  equation  can  be  used  to 

*  This  work  was  supported  by  the  Office  of  Standard  Reference  Data  of  the 
National  Bureau  of  Standards. 


2 

reproduce  all  the  available  P,p,T  experimental  data  as  well  as  to  produce 
limited  extrapolations  (based  on  thermodynamic  arguments)  of  the  surface 
into  important  regions  where  data  are  sparse.     The  range  of  the  equation 
is  bounded  at  low  temperatures  by  the  triple  point  temperature  (195. 48K) 
and  the  melting  curve  for  the  liquid,  and  at  high  temperatures  by  the 
isotherm  at  750K  (which  is  approximately  5/3  the  critical  temperature) . 
The  pressure  range  extends  to  8,000  bar.    No  mathematical  constraints 
were  imposed  on  the  equation,  and  only  P ,p ,T  data  were  used  in  the 
least  squares  fit  with  the  sole  exception  that  values  of  vapor  pressure 
were  explicitly  employed  in  regard  to  satisfying  conditions  of  phase 
equilibrium. 
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Following  Keenan  et  al.     the  Helmholtz  free  energy  function  was 
represented  in  reference  1  as  the  sum  of  two  terms,  the  first  being  the 
contribution  from  the  equation  of  state  and  the  second  a  function  of 
temperature  only  referring  to  the  properties  of  the  ideal  gas.     Thus  the 
Helmholtz  free  energy  was  expressed 

A(p,T)  =  A(p,T)  +  A°(T),  (1) 

where  A°(T)  is  the  contribution  of  the  ideal  gas.     T  is  the  absolute 
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temperature  in  kelvins  and  p  the  density  in  grams  per  cm  .    A  quantity 
Q(p,T)  was  defined  by 

A(p,T)   E  RT   [£n  p  +  pQ  (p,T)]  (2) 


Since         P  =  p     9A/dp , 
Eqs.    (1)  and  (2)  yield 

P  =  pRT  [1  +  pQ  +  p^  8Q/9p] . 


(3) 


Note  that 


Q  (p=0)  =  (3a) 


where  B2  is  the  second  vlrial  coefficient.     The  form  chosen  for  Q 


was 
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Q  =   ^      2]    a.  .  (r-ry-\  (4) 

1=1  j=l 


where  t  =  =  1.2333498  and  R  is  the  gas  constant.    Eq.   (3)  for 

the  pressure  was  fitted  In  the  least  squares  sense  to  the  experimental 
P,p,T  data.     The  results  of  this  fit  are  values  for  the  constants  a., 
listed  in  table  1.     (a^^  =  0  for  all  i,j  pairs  not  listed  in  table  1) . 
By  differentiation  of  eq  (1) ,  the  various  thermodynamic  quantities  can 
be  expressed  in  terms  of  Q: 
the  entropy. 


S(p,T)  =  R[£n  p  +  pQ  -  px  9Q/3x]  +  S°(T) ,  (5) 


the  internal  energy. 


E(p,T)  =  RpTx8Q/9x  +  E°(T),  (6) 


the  constant  volume  heat  capacity, 


C  (p,T)  =  -  Rpx^  9^Q/9x^  +  C  °, 


(7) 


Table  1 


-6.453022304053 
-13.719926770503 
-8.100620315713 
-4.880096421085 
-12.028775626818 
6.806345929616 
8.080094367688 
14.356920005615 
-45.052976699428 
-166.188998570498 
37.f 08950229818 
-40.730208333732 
1.032994880724 
55.843955809332 
492.016650817652 
1737.835999472605 
-30.874915263766 
71.483530416272 
-8.948264632008 
-169.777744139056 
-1236.532371671939 
-7812.161168316763 
1.779548269140 
-38.974610958503 
-66.922858820152 
-1.753943775320 
208.553371335493 
21348.946614397509 
247.341745995422 
299.983915547501 
4509.080578789798 
-37980.849881791548 
-306.557885430971 
24.116551098552 
-9323.356799989199 
42724.098530588371 
161.791003337459 
-507.478070464266 
8139.470397409345 
-27458.710626558130 
-27.821688793683 
298.812917313344 
-2772.597352058112 
7668.928677924520 


the  enthalpy  function. 


H(p,T)  -  RT(pQ  +  p2  aq/ap  +  px  3Q/3T  +  1]  +  E°  (T) ,  (8) 
the  constant  pressure  heat  capacity, 

Cp  <P»T)  -        +  R-  I  .  (8a) 

where, 

a  -  1  +  pQ  +  p2  3Q/3p  -    pT  3Q/aT  -  p^T  3^Q/3t3p  , 

and 

e-l+2pQ+4p2  3Q/3p  +  p32Q/3p2, 
the  heat  capacity  for  the  saturated  fluid  , 

^      (3P/3T)p  dPg 
"  S  "  ^    OP/3p)^    •    dT"  ' 

where       Is  the  vapor  pressure  of  the  liquid  • 

S**,  E**  and        are  the  corresponding  contributions  obtained  from  A°(T) 

A°(T)  -  (G°  -  E°)  -  RT(1  +  in  4.8180  T)  , 
S°--^A°(T)! 

,o.,o^,^.,dAl(T),  •  (10) 

-  -  T  d^  A°(T)/dT^, 


where  G°  is  an  analytic  representation  of  the  results  reported  by  Haar^  for 
the  properties  of  the  ideal  gas  state, 
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G  -  E 


RT 

i=2 


(lOa) 


where  is  the  energy  for  the  ideal  gas  at  0  K.  The  coefficients  a^  are 
listed  in  table  2. 

Table  2 


*1 

-3.872727 

.64463724 

^3 

3.2238759 

^ 

-.0021376925 

^5 

.86890833x10"^ 

*6  . 

-.24085149x10"^ 

*7 

.36893175xlO"-'-° 

^ 

^.35034664x10"-^^ 

a 

.2056303xl0"-'-^ 

mm 

-.6853420x10"^^ 

*11 

.9939243x10"^^ 

The  heat  capacity  values  and  the  other  thermodynamic  functions  calculated 
from  Eq. (10a)  for  the  temperature  range  100  K  <  T  <  1000  K  agree  with  those 
tabulated  in  reference  5  to  within  the  accuracy  of  those  values. 


Though    Eqs.   (1-lOa)  are  complete,  it  is  necessary  to  introduce  the 

Gibbs  phase  conditions  in  order  to  calculate  the  properties  for  the 

coexisting  phases.     However,  it  was  shown  by  Haar  and  Gallagher"'"  that 

almost  negligible  error  results  if  an  explicit  relation  is  used  for  the 

vapor  pressure,  P^,  obtained  by  separately  fitting  the  saturated  vapor 

pressure  data  of  Cragoe    for  the  range  from  the  triple  point  to  373K  and 
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of  the  mean  of  the  data  of  Beattie  and  Lawrence    and  of  Keyes    above  373K. 

The  equation  so  obtained  is 

P  T 

log^        =  ^    (AX  +  BX"^/^  +  CX^^^  +  DX^) 
c 

where         X  =  1  -  T/T 

c 

P    =  111.85  atm. 
c 

T    =  405. 4K 

c 

A  =  -7.296510 
B  =  1.618053 
C  =  -1.956546 
D  =  -2.114118 

Eq.   (10b)  and  Eq.   (3)  define  the  coexisting  phases  for  this  report 
and  all  properties  over  the  temperature-pressure  range  of  the  thermodynamic 
surface,  subject  to  the  low  temperature  boundary  of  the  melting  solid. 
The  thermodynamic  surface  is  consistent  with  the  following  values  for 
the  parameters  at  the  triple  point, 

T^  =  195. 48K 

P^  =  .06063  bar 

=  .00006382  g/cm^,        =  .73374  g/cm^, 

and  at  the  critical  point, 

T    =  405. 4K 
c 

P         113.04  bar 
c 

=  .2350  g/cm^ 


The  relationship  between  the  pressure  and  temperature  of  the  melting 

solid  was  calculated  by  means  of  the  Clapeyron  equation.     For  the  latent 
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heat,  the  value  reported  by  Overstreet  and  Giauque    was  used;  for  the 
specific  volume  for  the  solid  at  the  normal  melting  point,  the  value 
reported  by  McKelvey  and  Taylor''"^  was  used,  and  for  the  corresponding 
specific  volume  of  the  liquid  the  value  reported  by  Cragoe  and  Harper''"''" 
was  used. 

The  Clapeyron  equation  is  the  relation 


d|=ul^    dP.  (11) 


where  the  quantities  u'  and  u  are  the  specific  volumes  of  the  liquid 
and  solid,  respectively,  and  L  is  the  latent  heat  of  fusion.     From  the 
above  data,  the  quantity  — — r— ^    -  4x10  ^  atm  "'";  and  Eq.   (11)  can  be 
integrated  to  yield 


T  =  T    exp[4xl0  ^(P  -  P  )],  (12) 
s  s 


where  T    and  P    are  the  triple  point  values.     (The  differences  between  the 
s  s 

triple  point  values  and  those  of  the  normal  melting  point  are  negligible.) 
Also,  since  P^  at  the  triple  point  is  very  small,  the  relationship  can  be 
simplified  to 


T  =  195.48  exp  {4x10  ^  P  (atm)}. 


(13) 
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3.      Computer  programs. 

The  properties  calculated  from  the  equations  presented  in  this  report 
have  been  compared  with  the  various  thermodynamic  measurements  reported  for 
ammonia  in  the  comprehensive  review  in  reference  1.     In  almost  all  cases 
the  agreement  is  within  the  experimental  accuracy  of  the  data.  Comparisons 
of  results  with  P,p,T  data  are  given  here  in  figures  2  and  3.     It  was 
established  in  reference  1  that  for  most  of  the  vapor  phase  and  for  the 
coexisting  phases,  the  calculated  enthalpies  are  accurate  to  within  0.1%. 

To  facilitate  the  application  of  these  results  we  present  in 
Appendices  I  through  IV  computer  programs  in  FORTRAN  with  which  the  various 
properties  of  ammonia  can  be  calculated.     Four  such  programs  are  included: 

The  computer  program  in  Appendix  I  is  a  general  program  for  all  the 
thermodynamic  properties  discussed  in  section  2  of  this  report,  including 
the  program  for  the  calculation  of  Q  of  Eq.   (4)  and  its  derivatives  with 
respect  to  temperature  and  density  which  are  used  by  most  of  the  other 
programs.     The  independent  variables  are  either  pressure  and  temperature  or 
density  and  temperature. 

The  computer  programs  in  Appendix  II  refer  only  to  certain  properties 
for  the  coexisting  phases,  including  the  latent  heat  of  vaporization,  the 
vapor  pressure  of  the  saturated  liquid,  the  densities  of  the  saturated 
vapor  and  of  the  saturated  liquid  and  the  heat  capacity  of  the  saturated 
fluid.     The  dependent  variable  for  these  is  either  the  temperature  or  the 
corresponding  pressure  of  the  saturated  liquid. 
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As  all  of  the  above  programs  assume  densities  in  g/cm  ,  temperatures 
in  kelvins  and  pressure  in  atmospheres,  a  set  of  simple  routines  is  presented 
in  Appendix  III  for  conversions  back  and  forth  between  these  units  and  other 
commonly  used  systems  of  units. 
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Appendix  IV  contains  examples  of  typical  calculations  using  the 
routines  for  Appendices  I  through  III. 

All  quantities  other  than  pressure,  density  and  temperature  are 
calculated  and  used  in  dimensionless  units.     Results  in  two  commonly  used 
systems  of  units  can  be  obtained  as  follows: 

To  obtain  heat  capacities  and  entropy  in  SI  units  (joules/g),  multiply 
the  dimensionless  quantities  by  0.48820;  or  in  Btu/lb.,  multiply  the 
dimensionless  quantity  by  0.210027. 

To  obtain  enthalpies,  internal  energies  or  heats  of  vaporization  in 
SI  units  (joules/g  K) ,  multiply  the  dimensionless  quantities  by  (0.48820  T) 
where  T  is  in  K;  or  in  Btu/lb.  deg  F,  multiply  the  dimensionless  quantities 
by  (0.210027  T)  where  T  is  in  deg  R. 

The  programs  as  presented  are  in  single  precision  except  for  the 

coefficients  of  the  generating  function  'Q',  which  are  in  double  precision. 

This  will  give  accuracy  to  within  the  limits  of  the  equation  when  used 

on  most  computers.     (When  used  on  computers  with  word  lengths  of  48  bits 

or  more,  even  these  coefficients  may  be  used  in  single  precision.)     If  more 

precision  is  needed,  the  following  can  be  made  double  precision  in  all 

routines:  T,  P,  D,  QOl,  Q02,  QlO,  Q20,  Qll.     The  greater  precision  may 

be  required  when  making  calculations  at  liquid  densities  or  in  the  vicinity 

of  the  critical  point,  particularly  when  obtaining  values  for  the 

thermodynamic  functions  requiring  higher  derivatives  of  the  generating 

function  Q,  such  as,  for  instance,  C  .     In  these  cases  the  contributions 

p 

of  all  of  the  44  terms  are  of  importance,  as  there  is  much  cancellation 
of  contributions  from  individual  terms. 
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C  APPENDIX  I 

C  PACKAGE   OF   FORTRAN   PR3GRAMS   TO   GENERATS    THE  THERMO- 

C  DYNAMI    C   PROPERTIES  OF   AMMONIA  FOR   ENTIRE  REGION. 

C 

SUBROUTINE  QQ( Q«L • K» T* D I 

C      THIS    IS   THE   GENERATING        FUNCTION   FOR   Q   AND    ITS   DERIVATIVES  USED  IN 
C      THE  CALCULATION   OF   THE   OTHER        FUNCTIONS.      Q   IS   GENERATED   AS   A  FCT 
C     OF  T    AND   D,    AND  L    IS  THE  ORDER  DF   THE   DERIVATIVE    OF   Q   WITH  RESPECT 
C      TO  RECIPROCAL   TEMP,      AND  K   IS  THE   ORDER  THE   DSRIV.    W,R«T.  DENSITY. 

DOUBLE   PRECISION   A. Al , A2. DD. C . TT, QT 

DIMENSION   P(5,5)*0D( 10>,TT(7). AC 4A).II(AA} • JJf 44) • Al(  221* A2f 22) 

EQUIVALENCE   IA1U)>A<1)),(  A2(  1  )«  A(23)  ) 

DATA   I  1/6*1 .6*2 ,6*3 .6*4 ,4*5. 4*6, 4*7 •4*8» 4*9/ » N/44/ 

DATA  J J/1 •2*3«4«5,6,l«2«3.4,5,6»l»2*3*4«5*6.l.2«3»4,5*6*1.2«3»4«l« 
1    2*3.4*1*2«3*4«1«2«3«4» 1»2»3*4/ 

DATA  P/2«»2«»3..4.»5««2.«4.«6.»12..20«f»3»«6e«l2®e24^e«60,«4»«12.« 
1   24 ••48.«120. tS.«  20.»60.« 120..240./ 

DATA  Al /-6. 453022 3040500. -13. 719926T705D0, -8. 1005 2031  5700. -4 • 88009 
164  21D0*~12.0287756268D0.6.80634  593D0*  8. 08009436769D0. 1 4 • 356920006D 
20.-45.052976700*-166.18899857D0.37.90895023DO,-40.730  2083337DO, 

3  1.032994880700,55.84395581D0.492*016650817700.1737.9  36DO« 

4  -30.87491 5263766DO»71.48  353041600,-8.948264632DO,-169»7777441400« 

5  -1236.53237167D0.-7812.161 16831700/ 

DATA  A2/1 .7795482691 4D0. -3  8.9  7461 09600.-66.9  228588200,-1.75394377532 
13 200, 208. 553371 335O0.21348.9466D0 .247 .341 74600 .299. 98391 5550  0,4509 
2.  08057'379D0.-37980.  8500,-306.  5578854300.  24.  1 1 6551 100* -9323.  3568D0, 
34 2724. 0985306D0, I  61. 791 0033374600,-507. 478070464DO •9139. 4703974D0. 

4  -27453. 71 06300. -27^821 6888DD • 298.81 2917300,-2772.59735200, 

5  7668.92867800/ 

QsO  a 

IFCL+K)  12,14.18 
12  RETURN 

14  Ua500./T 

C=U- 1 . 233349778D0 

IF(0AaS{C) .LT.l.D-81  C=l.D-8 

TT<l)sl. 

DO   15  1=2,7 

15  TT<  I  )s!TT<  I-l  l*C 
lF(D.LT.l.E-8)    0*1. E-8 
DD(1)»1. 

DO   16  1=2.10 

16  DD( I )=DD<I-l )*D 
18  00   200  M=1.N 

I=II<M) 
J=JJ( M) 

IF(J-L-l)  2OO,2O»30 
20   IF(L)  200,24.28 
24  QT=1. 

GO   TO  50 
28  QT=P<L,L)/2. 

GO   TO  50 
30   IF<L)  200,32,40 
32  QT=TTCJ) 

GO  TO  50 
40  QT=P{L, J-l )*TTC J^LI 
50   IF«K)  200.60.70 


60  QT=QT*00( I > 

GO   TO  100 
70  Rsl. 

DO   72  MM=l,K 
72  R=R*I I-MM> 

QT=QT«R*OD(I-K) 
100  Q=Q«-A<M|4cQT 
200  CONTINUE 

RETURN 

END 

C 

SUBROUTINE  FZC T,FO, SO, CPO.UOl 
C     COMPUTES   THE   IDEAL   GAS        FUNCTIONS   AS  FCT  OF  T:   FREE   ENERGY! FO/RT ) 
C      ENTROPY( SO/R> ;      CPCCPO/R);  AND    INTERNAL   ENERGYC UO/RT I , 

DIMENSION  A(12) 

DATA   A/-3. 872727. .dA^ 63724. 3. 2 238759. 0021 376925, .8689 0833E- 5. 

1  240851 49E-7. .368931 75E-10,  -  . 35034664E-1 3  t . 20563027E-16 , 

2  -•6e5342E-20. •99392427E-24.0./ 

U=T 

IFCT.LE.IO.)  U*10. 

FO   =A(l  I^ALOGf  UI>A(2)/U4-A(  3l4-A(4l4cU-«-A(5l*U*U4-A<6)*U*«3^AC7)*U4t4t4 

1  4-A<8)«U**5*A«9>4«U**6*ACl  0  }*J**7*M  1  ll*U**8 

S0=A<1  )♦<  ALOGfUl^l.)    4-   A(3)    ♦   2.*AC4>*U   ♦  3**A(S}  *\J*\J   ♦  4.*A(6>*U 
1**3        5.*A(7>*U**4   ♦  6.*A(8)*U**5   ♦   7.*AC9>*U**6   ♦  8«*AC10)*U**7 

2  •»•  9.  *A<  11  )*U**8 

CVO=A<  1  )/U-l-2.*Af  4)4>6.*Af  5)*U4-12.*A(6)*U*U«-20.*AC7)*U**3-f30.*A(8)* 
1   U**4   ♦  42.*A(9)*U**5  ♦  56 . *Af 1 0 ) *U**6  ♦  72. *A< 1 1 ) *U**7 
CVO=-CVO*U 
SO=-S0 
U0=S0*F0-1 . 
RETURN 
END 

C 

FUNCTION  PRES<D.T> 
C      PRES   CALCULATES   THE  PRESSURE    IN  ATMDSPHERES   USING   O.T.Q   AND  DQ/D(D|. 
COMMON  /QQQQ/   QOO • 00 1 • QO 2 • Q 1 0 • Q20 « Q 1 1 
PRES   -   4.818*T*D»f l.*0*(Q00+D*Q01 ) ) 

RETURN 
END 

C 

FUNCTION  DPOD(D.T) 
C  OPDD  CALCULATES   THE  DERIVATIVE   OF   THE   PRESSURECIN  ATM)    WITH  RESPECT 
C  TO  DENSITY    CIN   GM/CC) . 

COMMON   /QQQO/  QOO • 001 • Q02 • 0 1 0 . Q20 . 01 1 

OPOO   =   4.818*T*< 1 •>2.*D*Q00   *  0*D*(4.*QOl^D*Q02} ) 

RETURN 

END 

C 

FUNCTION  ENTR<D.T.SO> 
C     CALCULATES  THE   ENTROPY    IN  01  MENS  I ONLESS   UNITS    (S/R).  REQUIRES 
C  PRIOR   CALL   TO  FZ,    QOO   AND  QIO 

COMMON  /QQQQ/  Q 00 . QO 1 . QO 2 • Q 1 0 • Q20 • Ql 1 

U=500 ./T 

TO=T*0 

lF<TD.LT.l.E-6)  TDsl.E-6 

ENTRsSIG't-D*U*Q10~O*Q00-ALOG(4.818*TO) 

RETURN 

END 


FUNCTION  ENTH(0«T.ENER) 
C     CALCULATES   THE  ENTHALPY    IN  DI MEMSIONLESS   UNITS    (H/RT).  REQUIRES 
C     PRIOR   CALL  TO  ENER*    FZ«Q00*Q01    AND  QlO. 

COMMON  /QQQQ/  000 • 001 • Q02* Q 1 0 • Q20 • 01 1 

ENTHsENER        1  •    ♦   ©♦QOO   ♦  D*D*QOl 

RETURN 

END 

C 

FUNCTION  ENER(D«T,UOI 
C     CALCULATES  THE    INTERNAL   ENERGY   IN   DI MENSIONLESS   UNITS  (U/RT>, 
C     REQUIRES  PRIOR  CALL  TO  FZ,    QOO   AND  QlO 

COMMON  /QQQQ/   QOO vQO 1« Q02* 01 0« Q20 • Ql 1 

U=500./T 

ENER=0*U*Q10   ♦  UO 

RETURN 

END 

C 

FUNCTION  CVCD«T,CPOI 
C     CALCULATES  THE  HEAT   CAPACITY   AT  CONST   VOLUME   IN  DI MENSIONLESS 
C      UNITS   (CV/RI*   REQUIRES  PRIOR   CALLS  TO  FZ.    QOO   AND  Q20* 

COMMON  /QQQQ/  QOO .QO 1 • Q02 • Ql 0 • Q20 fQl 1 

U=500./T 

C  V«CP  O-l • -0*U*U«Q20 

RETURN 

END 

C 

FUNCTION  CPCD.T.CV) 
C     CALCULATES  THE  HEAT  CAPACITY   AT  CONST  PRESSURE    IN  D I MENSI ONLESS 
C     UNITS   (CP/R).   REQUIRES  PRIOR  CALLS  TO  FZ«   QOO*    Q01«Q02«    QlO.  Q20* 
C      Qll    AND  CV* 

COMMON  /QQQQ/  QOO.QOl* Q02«Q10*Q20«01 1 

U=500./T 

CP={1  .+D*(  QOO-U4cQ10<»-D*(Q01-U*Ql  1  )  l>  ♦*2/<  1 .  ♦D*!  2«  *Q00€D*(  4.  ♦QO  l-^D* 
1  Q02)))^CV 
RETURN 
END 


SUBROUTINE  DFI ND( DOUT • P,D« T« OPDI 
C      THIS   SUBROUTINE,    BY   AN  ITERATIVE   PROCESS  USING   THE  -UNCTIONS 
C     PRES   AND   OPDD.    WILL  CALCULATE  THE   DENSITY  AS   A   FCT   OF   A  GIVEN 
C     PRESSURE    AND   TEMPERATURE.      THE   REQUIRED   INPUTS   A«IE   P.   T.  AND 
C     DC WHICH   IS  AN   INITIAL  GUESS  FOR  THE   DENSITYI.      THE  OUTPUTS  ARE 
C     DOUT   <THE  DENSITY  CDNSISTANT   WITH  THE    INPUTS  P   AMD  T>    AN9  DPD 
C      (THE   VALUE  OF   DP/DD  AT  P  AND   T>.      THE  VALUES  OF  P   ARE   IN  ATM. 
C     T   IN  K.   AND   D   AND  OOUT    IN   G/CC.      NOTE  THAT   FO^    T   LESS  THAN  TC. 
C      AND  FOR  P  AND  T   NEAR   THE  COEXISTANCE   CURVE.   CARE  MUST   BE   USED  IN 
C      THE   SELECTION  OF   THE   INITIA.    GUESS   F0«    O  SO   THAT  SPURIOUS  VALUES 
C      IN  THE   TWO  PHASE  REGION  ARE   NOT  RETURNED. 

COMMON  /QQQQ/   QO • Ql . Q2. QlO • Q20 • Ql 1 

DO=D 

L*0 
9  L«L4-1 

IF(OD.LE.0.)   00=1. E-8 

IFIDD.GT..85)  0D-.85 

CALL  QQ(QO.O.O.T.DDI 

CALL  QQCQl.O.l.T.OD) 

CALL  QQ(Q2.0.2.T.DD) 


DPDX=l • l*DPOO<DD«T) 
IFCOPOX.LT, .1  )  DPDX=,1 
PP=PRES«OD.T> 

IF< ABS( 1 .-PP/PI .LT.l .E-e)    GO   TO  20 
X={P-PP)/DPDX 

IF(DABS( X) .GT. • 1 )    X^X*. 1 /0A3SC X ) 

DD=OD-fX 

IF<L.Le.20l    GO  TO  9 
20  CONTINUE 
OOUT=00 
RETURN 
END 

C 

SUBROUTINE   PHASE( P»T. OL, DG.Kl 
C     THIS   SUBROUTINE   WILL   PERFORM   THE   FOLLOWING  FUNCTIONS: 


C  1.    IF   P   AND  T   ARE   SUPPLIED,    A   VALUE   OF   K   WILL   BE  RETURNED 
C  OF   -1,0«1«2  OR    3  WITH  THE  FOLLOWING   MEANINGS:    K-^l  IS 

C  ERROR   RETURN   ( P  OR   T  LESS  THAN   0.1        K=0   POINT   <P,T»    IS  ON 

C  COEXISTANCE   CURVE   AND   APPROXIMATE   VALUES   F0\    DL    AND   OG  ARE 

C  CALCULATED.   K=l    INDICATES  T  LESS  THAN  TC,   POINT    (P,T)  IS 

C  IN  VAPOR  PHASE.    APPROXIMATE   D=DG*OL  CALCULATED.      K«2  SIM- 

C  ILAR  TO   K=l.    BUT   POINT    « P»  T  I    IS    IN  LIQUID   PHIASE.  X«3 

C  INDICATES   T   G.T.    TC,    AND  D=DG=OL   IS  APPROXIMATED. 

C  2.    IF   T  LESS   THAN   TC   AND  P=0,    PS(T)    WILL   BE   CALCULATED  ALONG 
C  WITH  APPROXIMATE  DG   AND   DL .    IF        LESS   THAN  PC   AND  T=0, 

C  TS(PI    WILL    BE  CALCULATED   A.ONG  WITH  APPROXIMATE   DG   AND  DL. 

IFCP)  10.10.20 


10  IF(T.LT.195.    .OR.   T.GE.405.SI    GO   TO  90 

14  PsPSCT) 

16  CALL  DS(T«DL,OG| 

K  =  0 

RETURN 
20   IF(T)  22.22.30 
22   IF<P-lll.55l  24.24,90 
24    IFfP-.06)  90,26.26 
26  T=TS<P) 

GO  TO  16 
30   IF(T-405.5)  32.50,50 
32  PX=PS<TJ 
34  CALL  OS(T«DL,DGI 

IF(P-PX)  36,14,40 
36  DG=DG*P/PX 

DL«DG 

K=l 

RETURN 
40  DG=DL 
Ka2 

RETURN 

50  D=P«<  .00021  07-5.19O<-8*P-2.289D-7«TI-l.2823O>64'T 
Ds0-»-P/T/4.8l8 
IF<D.LT.1«D-SI  Ds.OOOOl 

DL=eD 
DG=OL 
K~3 
RETURN 
90  DL-O. 
OG^O. 

RETURN 
END 


c 

C  APPENDIX  II 

C  FORTRAN  PROGRAMS   FOR   FUNCTIONS  PERTAINING   TO  COEXISTANCE  LINE 

C 

FUNCTION  PS(T) 

C     PS  CALCULATES   THE   VAPOR   PRESSURE    (IN   ATM)   AS   A   FCT   OF   T   IN  K. 
C      FOR   T   GREATER   THAN   TC.    0.    IS  RETURNED. 

DIMENSION   A(4  )/-7» 296510,  1  .61  80  53* -1  •956546 « ->2 •  1  1  4  1 1  8/ 

IF(T.GT.405.4 )   60   TO  10 

TH=T/405«4 

TM=1 .-TH 

TS*0SQRT«TM) 

PsTM*TS 

Q=P4cTM 

RsQ4tQ 

F=<  A(l  )*TM+AC2>*P-I-A(3)*Q^AC  4»*R>/TH 
PS=11 1 .BS^EXPCF) 
RETURN 
10  PS=0. 
RETURN 
END 

C 

FUNCTION  OPSOTCT) 
C     CALCULATES  THE   DERIVATIVE   OF   THE    VAPOR   PRESSURE   CURVE  WITH 
C  RESPECT   TO  T.      FOR   T   GREATER   THAN   TC.    THE   VALUE   2.03    (THE  APPROX- 
C    IMATE   VALUE   AT   TC)    IS  RETURNED. 

DIMENSION  A(4)/-7,296510, 1 .dlSOSS*-! .956545 , -2, 1 1 4  1 1 3/ 

IF(T.GT. 405.4)    GO   TO  10 

TH=T/405.4 

TM=1.-TH 

TS=SQRT(TM> 

p*TM*TS 

Qsp*TM 

R«Q*Q 

F=(A(  I  )*TM'I-A(  2)*P+A(  3  )  *Q  +  A(  4)  ♦R  ) /TH 

PS=111.85«EXP(F) 

S=P*Q 

DPSDT5=(-F-A(  1  )-l.S*A(2)*TS-2.5*A(  3  )*P-5.«A (  4 )  *S )  ♦PS/T 
RETURN 
10  DPSDT=2.03 
RETURN 
END 

C 

FUNCTION  TSIP) 

C     CALCULATES  BY   AN    ITERATIVE   METHOD   USING   PSCT)    AND  GENERATING 
C      A   STARTING  VALUE   BY   AN  APPROXIMATE   EQUATION.    THE   SATURATION  TEMPERATURE 
C     AS  A   FUNCTION   OF   THE  PRESSURE    IH    ATM.      IF  P   LESS  TH^N   THE    TRIPLE  POINT 
C      VALUE   OR   GREATER   THAN  PC.    0    IS  RETURNED. 
Q=l ./P 

lF(P.GT..05e   .AND.   P.LT.111.85>    GO  TO  4 

TS=:0. 

RETURN 

4  CONTINUE 

T   =  252.10703  ♦  P*(5. 6658636  ♦  P*  C   -.12457077  «•   P*(. 00135987 
1    -P*. 00000531942 )) )    *   Q* C -1 5. 3 79678*Q* ( 1 .52831 5-0* .04906) ) 

5  PX*PS<T) 
IF(A8S(l.-P/PX).LT.l.E-5)    GO   TO  10 


OT  =  C P-PX l/OPSOTC T  » 
T»T-fOT 
GO  TO  5 
10  TS=T 
RETUPN 
END 

SUBROUTINE  OS(T»OL«OG} 
C      THIS   SUBROUTINE   WILL  CALCULATE    APPROXIMATE  COEXISTING  LIQUID 
C        AND   VAPOR  DENSITIES   AT  THE   TEMPERATURE   SUPPLIED,    USEFUL  F0(? 
C      INITIAL   GUESSES   IN   THE   SUBROUTINE   D=IN0«      IF  T    GREATER  THAN  TC, 
C      DL  =  DG  =  0. 

DC*. 235 

IF(T-405.4I  10,80,90 
10  DTsl .-T/405.4 

IF(T.LT.395» I    GO  TO  100 
T1*0T**#35 
T2«DT**»541 
T35:T2*T2 
T4=DT*2.9653 

DL=sT4*Tl*(  2«  117-1  •4097*T2-.89  802«T3I 
DG=T4-T1*( 2.1 17^1 •1390*T2*. 57253* T3) 
DL3DC*<DL*l .1 
DGsDC*(DG-H.I 
RETURN 
80  DL-OC 
DG=DC 
RETURN 
90  DL=0. 
DG=0. 
RETURN 
100  Q=DT 

DL=.3871  31->.00096947/Q4-Q*(  1  .1513 875 •»-Q*(-l  •4943106-l-Q*!*  118332  5)  ) 
DG=.  08288674-.00095e67/Q^Ql*<-.60039534-l-Q*<  1  .44  34594-0*1.  1678605)  ) 
RETURN 
END 

C 

FUNCTION  HV<T«DL.DG) 
C        THIS  FUNCTION  WILL  CALCULATE  THE   HEAT   OF  VAPORIZATION  IN  DIM- 
C        ENSIONLESS  UNITS   (L/RT)    AT   THE    INPUT   TEMPERATURE   T.      DL    AND  DG 
C        ARE   REQUIRED  TO   BE  PREVIOUSLY  CALCULATED  AND   SUPPLIED  AS  INPUT 
U=500./T 

CALL  QQ(OO.0.O.T,DL) 
CALL  QQf Ql «0*1 .T.DL) 
CALL  QQ(Q2*1 .O.T.OL) 
CALL  QQ(RO«0«0«T,DG) 
CALL  QQ<R1 .0. 1 .T. DG) 
CALL  QQ<R2. 1 .O.T.DG) 

HV=U*  (DG«R2-DL*Q2)    ♦   DG*(  R0»DG*R1  )    -   DL*<  Q0-I-DL*Q1  ) 

RETURN  ^ 
END 

C 

FUNCTION  CS(D.T.CPR) 
C     CALCULATES  THE   HEAT   CAPACITY   OF   THE   SATURATES   FLUID  IN 
C     OTMENSIONLESS   UNITS   (CS/R).    INPUTS  ARE   O.T   AND  THE  D I MENSI ONLESS 
C      CP/R   CALCULATED  PRIOR  TO  THE  CALL. 

COMMON   /QQQQ/   QOO  •  00  1  •  Q02  •  Q 1  0 •<320  •  Ql  1 

TDPDT=PRES(D.T) -0*500.  ^4.  81 8*0*  t  Ql  0-(>D*Ql  1  ) 

CSs=CPR-TOPDT/DPDD(D.T)/0/D   *DPSDT  (  T  )  * «  1 0  1 325 

RETURN 

END 


FUNCTION  TK(TF) 
TK=(TF+459,67)*5./9. 

RFTURN 
END 

FUNCTION  TF<TK) 
TF= 1 • 8*TK-A59.67 

RETURN 
END 

FUNCTION  PATM(PSIA) 
PATM=PSI A/1  4.69  6 
RETURN 
END 

FUNCTION  PSIA(PATM) 

PSIA=PATM*14«696 

RETURN 

END 

FUNCTION  PAT(P3AR) 
PAT=PBAR/1 ,01325 
RETURN 
END 

FUNCTION  PBAR(PATM) 
PB>6R=P  ATM*  1.0  1325 
RETURN 
END 

FUNCTION  DGCC(DLRCF) 
DGCC=DLBCF/6  2 .428 
RETURN 
END 

FUNCTION  OLBCF(DGCC) 
DLBCF=DGCC*62 .4  28 
RETURN 
END 


APPENDIX  IV 

EXAMPLES   OF   USE   OF   PROGRAMS   IN  ABOVE  PACKAGES 

THE  F0LL3WING   SAMPLE   PROGRAVI    WILL   CALCULATE  THE   DENSITY   AND   FIVE  THERMO- 
DYNAMIC  FUNCTIONS    AT   A   POINT   <P,T)    SUPPLIED,      P   SUPPLIED    IN   BAR,   T    IN   DEG  C, 
DENSITY   IN  GM/CC,   CPtCV.S   IN   JOULES/GM   DEG   C,      AND   INTERNAL   ENERGY,  ENTHALP") 
IN  JOULES/GM. 

COMMON  /QQQQ/   QOO • 00 1 • Q02 « Q 1 0 • Q20 • Ql  1 

1  READ   2,  PBAR.TC 

2  FORMAT<2F10,3,Fl0.6,2F10.3,3Fi0.5) 

PsPATM(PBAR| 
IF<P,LE.O,J  STOP 

T=TC^273,15  * 

CALL  PHASE<P,T,OL,OG«Kl 

IF(K.EQ.-1 )    GO  TO  I 

CALL   DFIND( D,P,DL,T,OX) 

CALL   QQ(Q10,1 .O.T.DI 

CALL   QOf 020,2,0,T,D) 

CALL   0Q(Q11,1 ,1 .T,0) 

CALL  FZ(T,FO,SO,CVO,EO) 

CVV=CV(D,T.CVO) 

CPP=CP(D,T,CVV)*.4882 

CVV=CVV*.4882 

S=ENTR(0,T,SO)* .4882 

E=eNER(0,T,EO) 

H=ENTH(D,T,E) *.4882*T 

E=E*.4882*T 

PRINT  2,PBAR,TC,D,H,E,S,CPP,CVV 

GO   TO  I 
END 


A    SAMPLE   PROGRAM  TO   CALCULATE   THE   SATURATION  TEMPERATURE 
IN   DEG  F,      DENSITIES    IN   LB/CU   FT,    AND   HEAT   3F   VAPORIZATION  IN 
BTU/LB  AT   A   GIVEN   PRESSURE    IN  PSIA 

READ   1.  PSI 
1  FORMAT(5F10.3> 
PsPATM(PSI ) 
T=TS(P) 

CALL   OSCT, XDL.XDG) 

CALL  DFlND(DL,P.XOL,T,OX) 

CALL  DFIND(DG,P,XDG,T,DY) 

HH=HV(T,OL,DG»«.210027*T 

OL=OLBCF(OL> 

0G=DL8CF(DG) 

PRINT    1,    PSI ,TT.OL,OG,HH 

STOP 
END 
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